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The  exact  solution  for  the  dynamic  stress  and  strain  waves  generated  by  a 
point  explosive  source  in  a  whole  space  is  a  simple  linear  combination  of  the 
reduced  displacement  potential,  the  displacement  pulse  and  the  velocity 
pulse.  The  first  decays  as  R'3,  the  second  as  R~2,  and  the  third  as  R*1 
(Haskell  1964).  The  far- field  strain  wave  pulse  can  be  generated  from  the 
velocity  pulse  by  simply  multiplying  the  latter  by  the  slowness  of  the  medium. 
This  suggests  that  a  similar  close  relationship  between  the  dynamic  strain  and 
stress  fields  and  the  dynamic  velocity  field  might  exist  in  more  general 
media.  Large  data  bases  of  near- field  velocity  recordings  have  been  collected 
and  analyzed  in  the  past  whereas  high  frequency  dynamic  strain  measurements 
have  seldom  if  ever  been  made.  The  purpose  of  this  study  has  been  to  attempt 
to  use  the  data  base  of  velocity  records  to  infer  the  levels  of  dynamic  strain 
and  stress  that  typically  occur  near  to  explosions  and  earthquakes. 

The  reason  that  near -field  dynamic  strains  are  of  current  interest  is 
that  a  variety  of  laboratory  studies  have  shown  that  some  typical  geological 
materials  begin  to  behave  in  a  nonlinear  fashion  at  shear  strain  levels 
between  10s  and  10'*  (Havko,  1979).  This  nonlinearity  manifests  itself  as  a 
lowering  of  effective  Q.  The  phenomenon  is  frequency,  amplitude  and 
overburden  pressure  dependent.  Day  and  Minster  (1986)  have  discussed  the 
potential  significance  of  nonlinear  behavior  on  interpretation  of  ground 
motion  data  from  a  series  of  small  explosions  in  salt,  but  the  events  were  so 
small  that  It  Is  not  straightforward  to  relate  their  results  to  nuclear 
explosions  or  earthquakes.  Also,  the  recording  geometry  of  their  experiment 
was  so  simple  that  they  could  use  whole  space  wave  propagation  results. 

Near -field  wave  propagation  is  usually  much  more  complex  near  to  general 
seismic  sources  (Burdick  et.  al . ,  1986,  Heaton  and  Helmberger,  1978)  It  is 
important  to  establish  the  levels  of  dynamic  strain  that  occur  for  these  more 
typical  sources  in  media  where  wave  propagation  is  more  complex 


We  present  here  the  results  of  an  Investigation  of  the  relationship 
between  the  velocity  and  strain  fields  near  to  seismic  sources  in  realistic 
crustal  structures.  First  we  develop  the  theory  necessary  for  computing 
stresses  and  strains  (rather  than  seismic  motions)  using  the  generalized  ray 
methodology.  Then  we  demonstrate  that  in  a  layered  half  space  the  far  field 
strain  wave  is  indeed  very  similar  to  the  far  field  velocity  wave.  We  develop 
a  transfer  operator  which,  when  convolved  with  an  observed  velocity  record, 
produces  an  estimate  of  a  strain  record.  As  noted  above,  in  a  whole  space, 
this  operator  would  be  a  delta  function  with  the  amplitude  of  the  slowness  of 
the  medium.  Ue  show  that  in  a  layered  half  space  it  is  a  delta- like  function 
with  approximately  the  amplitude  of  the  slowness  at  the  receiver.  (For  some 
strain  components,  it  is  the  amplitude  of  the  vertical  slowness  and  for  others 
the  horizontal  slowness.)  Next  we  apply  these  transfer  functions  to  a  suite 
of  velocity  records  from  nuclear  explosions  and  from  a  small  earthquake  to 
establish  the  levels  of  dynamic  stress  and  strain  which  they  generated.  The 
strain  levels  for  the  earthquake  were  found  to  be  comparable  to  those  for  a 
modest  nuclear  explosion  once  the  difference  in  source  depths  is  accounted 
for.  Presuming  that  the  laboratory  measurements  of  nonlinearity  are  relevant 
to  the  earth,  the  results  of  this  study  indicate  that  there  is  a  thin  region 
at  the  surface  of  the  earth  where  effective  Q  can  be  expected  to  be  very  low. 

STRESSES  &££  STRAINS  Hi  A  WHOLE  SPACE 

The  strain  tensor  in  a  whole  space  for  a  point  explosive  source  has  only 
one  nonzero,  far  field  component  in  spherical  coordinates.  This,  the  radial 
component,  is  given  by 

F  dV  „/  dR~  -2p/ R3  -2v' /aR*  -  p" /a2  R  (1) 

The  radial  stress  is  given  by 

P  -  -  lpp3{p/  R3  -  p- /aR2  ]  -  pv  /  R  O) 


•  is  the  reduced  displacement  potential,  a,  0  and  p  are  the  compressional 
velocity,  shear  velocity  and  density  of  the  medium  and  the  primes  denote 
differentiation.  The  last  term  on  the  right  in  Equation  1  is  just  the  far 
field  velocity  pulse  divided  by  the  compressional  velocity.  Thus,  the  strain 
and  stress  fields  are  given  exactly  in  terms  of  linear  combinations  of  the  RDP 
and  its  derivatives.  Seismic  motion  parameters  generally  have  a  near-field 
term  which  decays  as  R'2  and  a  far-field  term  decaying  as  R*1-  Strains  and 
stresses  have  additional  near  field  terms  that  decay  as  R"3. 

Of  course,  these  expressions  become  much  more  complicated  for  a  general 
point  source  in  a  layered  half  space  where ,  among  other  complications , 
cylindrical  coordinates  must  be  used.  However,  these  half  space  expressions 
become  much  more  tractable  if  near- field  terms  can  be  neglected.  The  simple 
whole  space  expressions  can  be  used  to  estimate  over  what  ranges  the 
near-field  terms  are  significant.  Figure  1  shows  the  radial  strains  for  two 
point  explosions  with  yields  of  100  and  1000  kt.  The  RDP’s  were  computed 
using  the  formalism  of  Helmberger  and  Hadley  (1981)  with  the  yield  scaling 
relations  published  by  Burger  et  al.  (1987).  The  medium  parameters  are  given 
in  Table  1.  The  far-field  term  is  shown  in  the  top  row,  the  near-field  terms 
in  the  second  row  and  the  sum  at  the  bottom.  At  a  range  of  1  km,  the 
near- field  terms  are  important.  By  3  km  however,  they  have  become  relatively 
insignificant.  This  result  appears  to  have  little  dependence  on  yield.  In 
the  next  section,  we  will  present  a  theory  for  computing  body  wave  stresses 
and  strains  using  only  the  far-field  terms.  Based  on  this  calculation,  it 
should  be  valid  for  ranges  larger  than  3  km. 

The  amplitudes  of  the  strain  pulses  shown  in  Figure  1  should  give  a 
preliminary  indication  of  how  significant  the  nonlinear  material  behavior 
observed  by  Mavko  (1979)  is  likely  to  be.  Figure  2  has  been  redrafted  from 
his  paper  and  illustrates  his  basic  result.  The  value  of  Q  is  observed  to 
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Figure  Lj.  The  analytic  solution  for  the  radial  strain  generated  by  a  point 
explosion  In  a  whole  space.  The  left  column  Is  for  a  range  of  1  km.  and  the 
right  for  a  range  of  3  km.  An  RDP  for  a  100  kt.  explosion  is  shown  at  the  top 
and  for  a  1000  kt.  shot  at  the  bottom.  In  each  Instance  the  near  field  terms 
are  shown  first  by  themselves,  then  the  far  field  term  and  then  the  sum  or 
complete  solution.  The  near  field  terms  become  insignificant  by  3  km. 
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begin  to  increase  at  about  10*6  strain.  By  10*5  strain  the  effect  is  very 
significant  in  most  cases.  The  magnitude  of  the  effect  appears  to  also  depend 
strongly  on  the  geologic  material  involved.  It  is  generally  believed  that  the 
decrease  in  Q  is  caused  by  frictional  sliding  on  cracks  (Stewart  et  al., 

1983).  The  observed  effect  should  therefore  also  depend  on  the  density  of 
cracks  in  the  material.  As  overburden  pressure  closes  these  cracks,  it 
becomes  increasingly  difficult  for  sliding  to  occur.  Thus  the  nonlinear  zone 
probably  does  not  extend  to  any  great  depth  in  the  earth.  Extending  the  1  km 
result  in  Figure  1  out  in  range  by  dividing  by  R,  we  find  that  for  a  1000  kt 
explosion,  far-field  strain  drops  to  10'5  at  about  60  km  and  does  not  fall  to 
10'6  until  almost  600  km  The  nonlinear  zone  would  appear  to  occur  over  a  very 
significant  region. 

STRESSES  AND  STRAINS  IN  £  LAYERED  HALF  SPACE 

Once  a  far-field  approximation  has  been  made  in  generalized  ray  theory, 
the  differences  in  the  expressions  for  seismic  motions  and  seismic  strains  are 
actually  quite  small.  Stress  tensors  can  be  generated  from  strain  tensors 
through  the  usual  definitions.  The  necessary  modifications  can  be  most  easily 
demonstrated  by  beginning  with  the  expressions  presented  by  Langston  and 
Helmberger  (1975).  Figure  3  shows  the  coordinate  system  in  which  their 
formalism  is  based.  The  small  rectangle  at  the  origin  represents  a  fault  with 
given  dip  (6)  and  rake  (A).  Here  we  will  consider  the  possibility  of  explosive 
sources  as  well.  The  vertical,  radial  and  tangential  motions  are  given  by  W, 

Q  and  V  respectively.  The  displacements  are  given  in  terms  of  the  Laplace 
transformed  seismic  potentials  n  and  X  by  the  Langston  and  Helmberger 
(1975)  equation  5 


where  the  near -field  terms  they  give  have  been  dropped.  The  Laplace 
transformed  potentials  for  a  generalized  ray  are  given  in  their  equation  6  by 


(4) 


A/0 

4/rp 


M0  is  the  seismic  moment,  p  is  the  density  at  the  source  and  p  and  r)  are  the 
horizontal  and  vertical  slownesses  respectively.  The  phase  of  the  ray,  P,  is 
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where  the  thj's  are  the  thicknesses  of  the  layers  encountered  by  a  given  ray. 
The  AVs  are  horizontal  radiation  pattern  terms  and  the  Cj's,  SVj's  and  SHj's 
are  vertical  radiation  pattern  terms.  The  J  indices,  0  through  3,  denote  an 
explosion,  a  vertical  strike  slip,  a  vertical  dip-slip  and  a  45  °  dipping 
normal  fault  viewed  at  an  azimuth  of  45*  respectively.  Any  other  point  source 
can  be  built  from  these  fundamental  sources  through  linear  combination.  These 
radiation  terms  are  not  of  importance  to  this  discussion,  so  they  will  not  be 
explicitly  written  out  here.  They  are  given  in  Langston  and  Helmberger 
(1975),  though  we  have  added  the  index  0  for  an  explosion  source  which  they 
did  not.  For  an  explosion,  A0  is  1  and  C0  is  1/a2.  The  factor  FI  represents 


the  product  of  transmission  and  reflection  coefficients  along  the  raypath. 

The  inverse  transform  is  evaluated  using  the  usual  Cagniard  de-Hoop  method 
(Wiggins  and  Helmberger,  1974,  Helmberger,  1974).  It  is  important  to  follow 
in  some  detail  how  the  displacements  in  Equation  3  are  related  to  the 
potentials  given  in  Equation  4.  The  spatial  derivatives  in  Equation  3  can  be 
evaluated  prior  to  the  inverse  transforms.  If  near  field  terms  are  neglected, 
they  simply  produce  some  additional  multiplicative  factors  in  Equation  4. 
Helmberger  (1974)  discussed  these  factors  in  some  detail,  and  we  shall  follow 
his  notation  and  utilize  many  of  his  results  here.  We  need  to  derive  these 
factors  for  both  the  case  in  which  the  receiver  is  embedded  in  the  half  space 
and  when  it  is  on  the  free  surface.  We  shall  show  that  the  only  difference 
between  seismic  motion  generalized  rays  and  strain  generalized  rays  is  that 
they  have  slightly  different  multiplicative  factors,  hereafter  called  receiver 
functions . 

To  illustrate  how  these  receiver  functions  are  derived,  let  us  consider 
the  case  shown  in  Figure  4.  We  have  an  incident  P  generalized  ray  and 
receiver  points  buried  at  depth  h  in  the  top  layer  of  the  stratified  half 
space.  The  arriving  P  ray  carries  with  it  a  4>  potential  with  initial 
amplitude  A0  If  we  wish  to  compute  the  vertical  component  of  motion,  we 
substitute  the  expression  for  the  ^  from  Equation  4  into  the  expressio:  for  W 
In  Equation  3.  For  a  single  ray  In  the  ray  sum,  the  result  is  of  the  foi 

«• 

W- j  sR„A0dp  (6) 

o 

A'  BZ  m~€'rja 

The  factor  of  s  generates  a  time  deriv.  ’-.ive,  and  Rpz  is  the  receiver  function. 
The  first  of  the  subscripts  indicates  tha  the  initial  potential  was  for  a  P 
wave.  The  second  indicates  that  the  vertical  component  of  motion  was 
computed.  The  converted  wave  shown  at  the  top  c  Figure  4  would  have  produced 
a  similar  result  except  that  the  receiver  function  -ould  be  R^.  Also,  the 
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product  n  would  include  the  free  surface  conversion  coefficient  Rn 
reflected  P,  the  receiver  function  would  again  be  Rpt  but  the  value  of  the 
function  would  change  sign  due  to  the  factor  t'  .  It  is  defined  to  be  -1  if 
the  last  leg  of  the  ray  is  upgoing  (Incident  P)  and  +1  if  the  last  leg  of  the 
ray  is  downgoing  (reflected  P) .  It  arises  because  the  phase  of  the  ray  is 
increasing  with  z  in  one  case  and  decreasing  in  the  other.  It  is  prised  to 
distinguish  it  froa  the  siailar  tens,  «  defined  by  Langston  and  Helaberger 
(1975)  which  has  siailar  properties  depending  on  whether  or  not  the  ray  is  up 
or  downgoing  at  the  source  rather  than  the  receiver .  As  indicated  at  the  top 
of  Figure  4,  the  free  surface  reflection  coefficient,  would  again  be 
included.  The  coaplete  set  of  receiver  functions  for  receivers  buried  in  a 
layer  are  given  at  the  top  left  of  Table  2.  It  shows  that  converting 
potentials  to  dlsplaceaents  involves  no  aore  than  aultiplylng  by  plus  or  minus 
the  vertical  or  horizontal  slowness  of  e  ray 

Ue  next  consider  the  coaputatlon  of  strains  within  the  aediun  using 
generalized  ray  theory  To  transfora  the  three  coaponent  displaceaent ,  D, 
into  a  strain  we  need  to  take  spatial  derivatives  according  to  the  fora 

£„  -  ^dD,/dxt*  dD,/dxt  '  (  /  ) 

Ue  can  again  take  these  derivatives  of  Equation  3  prior  to  the  inverse  Laplace 
transfora  Again,  they  will  result  in  additional  aultiplylng  factors  of  sp 
or  -s«n  Just  as  in  the  whole  space  case,  transforaing  from  seismic  motion 
to  strain  requires  only  aultiplication  by  slowness  The  additional  facte:  ct 
s  causes  a  time  derivative  which  transforms  the  displacement  components  to 
velocity  The  strain  and  velocity  tiae  histories  appear  to  be  closely  related 
just  as  in  the  whole  space  case  The  additional  slowness  aultipliers  will 
cause  phase  shifts  in  the  Inverse  Laplace  transforas  so  the  pulses  are  no 

longer  exactly  proportional  Ue  shall  show,  however,  that  in  most  cases  the 
effect  of  the  phase  shifting  is  not  large  The  strain 
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receiver  functions  for  receiver  points  within  the  aediun  ere  given  at  the 
bottoa  left  of  Table  2.  A  natural  extension  of  the  notation  of  Helmberger 
(1974)  has  been  used.  The  contribution  to  the  zz  component  of  the  strain 
tensor  due  to  an  incident  P  wave  is  computed  using  the  coefficient  Rpzz  For 
an  incident  S  wave,  it  would  be  RgU  and  so  on.  There  are  10  strain  receiver 
coefficients  and  only  5  motion  coefficients,  but  as  shown  in  Table  2  they  are 
very  closely  related  and  they  appear  in  Equation  6  in  an  identical  fashion 

Most  velocity  records  froa  explosions  and  earthquakes  come  from 
instruments  located  on  the  surface  of  the  earth.  It  would  thus  be  best  if  we 
could  establish  something  about  the  levels  of  stress  and  strain  in  the 
near-field  utilizing  this  information.  Furthermore,  if  we  can  develop  models 
that  successfuly  predict  levels  of  strain  and  strain  decay  rates  at  the 
surface  of  the  earth,  we  can  then  use  these  models  with  some  confidence  to 
predict  strain  levels  and  decay  rates  within  the  earth.  To  derive  the  seismic 
motion  receiver  functions,  we  determine  the  composite  response  of  the  three 
phases  illustrated  at  the  top  of  Figure  4  in  the  limit  that  h  becomes  small 
compared  to  the  shortest  wavelength  of  interest.  These  are  essentially  factors 
which  are  the  sums  of  the  amplitudes  of  the  direct  and  reflected  waves  near 
the  free  surface  The  geometry  for  the  incident  P  wave  case  is  illustrated  in 
Figure  4  When  the  P,  pP  and  pS  arrive  coincident  in  time  the  composite 
amplitude  of  the  response  is  just  A0(RPZ+RppRpz+RpsRsz)  The  reflected  waves  for 
incident  SV  would  be  SS+SP  and  for  SH  Just  SS .  This  means  there  is  some  finite 
depth  under  a  surface  velocity  receiver  (defined  by  the  shortest  wavelength  oi 
interest  criterion)  where  the  strains  estimated  from  that  receiver  can  be 
thought  of  as  existing  The  highest  frequency  of  interest  here  is  about  hr 
and  the  slowest  wave  speed  about  1  km/sec  Thus,  the  strains  that  we  estimate 
probably  exist  in  the  earth  down  to  depths  of  200  m  or  more 

The  final  step  needed  is  to  derive  the  receiver  coefficients  for  strain 
at  the  free  surface  It  may  appear  that  a  reasonable  wav  to  accomplish  this 


U 


would  be  to  dif ferentlace  Equation  6  while  inserting  the  spproprlete 
displacement  receiver  coefficient  end  to  make  the  usual  far  field 
approximation .  However,  this  approach  leads  to  nonzero  vertical  stresses  on 
the  free  surface.  To  obtain  the  correct  result,  it  is  important  to  first  take 
the  required  spatial  derivatives  of  the  displacements  generated  by  the  three 
phases  shown  at  the  bottom  of  Figure  4  and  then  to  take  the  limit  as  h  becomes 
small.  The  changes  in  sign  caused  by  the  factor  c'  and  the  factors  of  t)e 
generated  by  vertical  derivatives  of  the  S  wave  potential  instead  of  qa 
generated  by  the  P  wave  potential  cause  a  very  significant  difference.  The 
correct  receiver  functions  for  computation  of  free  surface  strains  are  shown 
at  the  bottom  right  of  Table  2.  They  are  factors  of  •  slowness  and  •  a 
dimensionless  constant  times  different  seismic  motion  receiver  coefficients. 
Note,  however,  that  for  motions  in  the  medium  Rrzz  is  generated  from  Rpz 
whereas  at  the  free  surface  it  is  generated  from  R^  Several  other  similar 
changes  exist.  The  cancellation  of  the  terms  Rrzx  and  Rraz  along  with  Rj^  and 
RsrZ  guarantees  that  Ea  will  be  zero.  The  OU/30)  has  no  far  field 
contribution,  and  ^  is  identically  zero,  so  Eu  will  be  zero.  The  diagonal 
term  Ezz  is  a  linear  combination  of  the  terms  Rm,  Rm,  R^  and  R^  which 
always  remains  zero  though  the  algebraic  details  are  somewhat  more 
complicated.  The  free  surface  condition  is  satisfied  by  the  strain  receiver 
coefficients  at  the  bottom  right  of  Table  2. 

The  preceding  has  shown  that  the  relationship  between  velocity  and  strain 
in  a  half  space  is  very  similar  to  their  relationship  in  a  whole  space. 

Instead  of  multiplying  velocity  by  the  composite  slowness  of  the  medium  to 
obtain  strain,  however.it  is  necessary  to  multiply  by  either  the  vertical  or 
horizontal  slowness.  Figure  5  illustrates  the  degree  to  which  this  is  true 
for  a  homogeneous  half  space.  The  source  is  the  same  1000  kt  explosion  used 
in  Figure  1,  and  tne  medium  parameters  are  again  given  in  Table  1  On  the 
left  are  the  velocity  traces  Each  generalized  ray  in  the  sum  has  been 
multiplied  by  a  constant  value  of  slowness  equal  to  the  value  of  p  or  »?  at 


lAtOLt  A  coapar i  son  of  velocity  trace*  multiplied  bv  constant  values  of 
slowness  (left)  with  actual  strain  traces  (right)  The  traces  compare  well 
when  the  receiver  is  buried  Differences  occur  when  the  receiver  is  on  the 
free  surface  due  to  the  free  surface  condition 
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which  the  Cagniard  contour  for  that  ray  leaves  the  real  axis.  On  the  right 
are  the  exact  strains.  The  case  of  a  surface  receiver  is  shown  on  the  top  and 
of  a  receiver  buried  at  15  km  on  the  bottom.  In  each  case,  the  range  is  15 
km.  At  the  surface,  the  shape  of  the  velocity  times  slowness  pulse  shape  is 
virtually  the  same  as  the  Ezz  shape.  The  difference  in  amplitude  is  due  to 
the  fact  that  Rpzz  is  not  just  fja  Rpz  but  the  more  complicated  form  given  in 
Table  2.  If  the  velocity  pulse  is  multiplied  by  the  correct  factor,  again 
evaluated  at  the  constant  value  of  p,  then  the  amplitudes  of  the  pulses  agree. 
The  pulses  agree  in  shape  and  amplitude.  The  true  Ej^  strain  is 
identically  zero  which  is  not  predicted  by  the  approximation.  Within  the 
medium,  the  approximation  works  very  well.  The  reason  it  does  not  work 
exactly  is  because  of  the  phase  shifting  that  occurs  as  p  evolves  along  the 
complex  Cagniard  contour. 

Though  the  effect  of  the  phase  shifting  in  Figure  5  is  small,  a 
homogeneous  half  space  is  a  very  simplified  medium.  In  more  realistic  crustal 
structures  where  unusual  raypaths  may  be  important,  the  phase  shifts  can  have 
a  large  effect.  It  is  necessary  to  find  a  way  to  account  for  this  phase 
shifting  no  matter  whether  it  is  strong  or  weak.  The  most  straightforward 
approach  to  devising  such  a  velocity  to  strain  transfer  method  is  to  utilize 
theoretical  frequency  dependent  transfer  operators.  These  operators  are 
generated  by  computing  a  theoretical  velocity  and  a  theoretical  strain 
response  for  a  site  and  deconvolving  the  former  from  the  latter.  Estimates  of 
strain  records  are  generated  from  velocity  records  by  convolving  the  velocity 
records  with  the  transfer  operators.  We  shall  illustrate  in  the  following 
that  in  most  cases  these  transfer  operators  are  very  delta-like  functions  with 
amplitudes  controlled  by  the  velocity  at  the  receiver  site.  The  latter  is  a 
relatively  well  known  quantity  and  we  believe  that  as  long  as  the  frequency 
shift  between  the  velocity  and  the  strain  records  is  small  the  transformed 
records  should  be  good  estimates  of  actual  strain  time  histories 


STRESSES  AND  STRAINS  NEAR  IQ  NUCLEAR  EXPLOSIONS 


The  first  data  base  to  be  processed  using  the  transfer  operators  is  a 
suite  of  near  field  velocity  recordings  from  five  NTS  nuclear  explosions.  The 
events,  which  were  detonated  at  Pahute  Mesa,  are  SCOTCH  (155  kt),  INLET  (324 
kt),  MAST  (406  kt) ,  ALMENDRO  (670  kt)  and  BOXCAR  (1300  kt) .  (The  yields  are 
taken  from  Burger  et.  al.  (1987).)  The  recording  sites  were  at  horizontal 
ranges  between  3.3  and  22.5  km.  The  recording  instruments  were  L7  velocity 
meters  which  have  a  response  flat  to  velocity  throughout  the  seismic  band,  so 
the  signals  recorded  on  them  are  essentially  velocity  versus  time.  The  event 
locations,  recording  lines  and  receiver  locations  are  shown  on  a  map  of  Pahute 
Mesa  in  Figure  6.  Not  all  of  the  available  records  were  suitable  for  use  here 
since  some  were  too  close  and  since  some  recording  channels  failed. 

Each  recording  site  used  needed  to  have  produced  a  three  component  data 
set.  In  order  to  compute  the  transfer  functions,  it  is  also  necessary  to  have 
theoretical  estimates  of  the  velocity  and  strain  pulses  for  each 
source -station  pair.  The  methodology  for  computing  such  synthetics  is 
discussed  in  Helmberger  and  Hadley  (1981)  and  in  Burdick  et  al .  (1984). 
Computation  of  body  waves  which  are  the  phases  of  interest  here  is 
accomplished  through  summation  of  generalized  rays.  Several  possible  plane 
layered  models  for  the  crustal  structure  at  Pahute  Mesa  are  available  in  the 
literature,  but  not  all  of  them  produce  synthetic  near  field  seismograms  which 
closely  match  the  observations.  The  one  used  in  this  study  was  provided  by  S. 
H.  Hartzell  (personal  communication)  and  was  specifically  developed  to  produce 
accurate  near  field  synthetics.  It  is  shown  in  Figure  7  along  with  the 
alternate  models  of  Helmberger  and  Hadley  (1981),  Hamilton  and  Healy  (1969) 
and  Carroll  (1966).  The  three  models  are  in  basic  agreement,  differing  only 
in  the  fine  detail  of  the  gradients.  The  parameters  of  the  model  are  given  in 
Table  1 . 
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Figure  7 .  p  and  S-wave  velocity  structure  for  Pahute  Mesa  obtained  in  this  study  (solid 
curve)  compared  with  two  other  proposed  structures  for  the  area  (dashed  and  dotted  curves) 


An  example  of  the  computation  of  a  transfer  function  for  the  Pahute  Mesa 
crustal  structure  is  given  in  Figure  8.  The  calculation  is  for  the  same 
theoretical  1000  kt  event  considered  in  Figure  1  and  Figure  5.  It  was  placed 
at  a  depth  of  1  km  in  the  crust.  The  synthetic  vertical  and  radial  veloc  ity 
pulses  are  shown  on  the  left,  the  four  nonzero  partial  derivatives  of  velocity 
with  respect  to  spatial  coordinates  in  the  center  and  the  transfer  functions 
are  on  the  right.  The  generalized  ray  sum  used  in  computing  the  medium 

response  contained  the  primary  rays  of  the  P  and  pP  type  and  near  receiver 

conversions  of  the  PS  type.  The  amplitude  of  the  P  wave  velocity  pulse  is 

about  10  cm/s,  the  strain  amplitudes  are  a  few  tens  of  ^strain  and  the 

transfer  functions  a  few  hundreds  of  pstrain/cm.  The  top  strain  is  Ezz  and 
the  bottom  is  The  center  two  must  be  summed  to  form  Note  that  they 

sum  exactly  to  zero  which  guarantees  the  conservation  of  the  free  surface 
condition.  The  transfer  functions  are  computed  by  transforming  the  velocity 
pulse  and  corresponding  partial  derivative  pulse  into  the  frequency  domain 
using  a  fast  fourier  transform  algorithm.  The  latter  is  divided  by  the  former 
and  the  inverse  transform  taken.  The  resulting  transfer  operator  strongly 
resembles  a  delta  function  with  a  signal  to  noise  ratio  of  better  than  5  to  1 . 
To  suppress  the  noise  further  a  gaussian  filter  with  a  cutoff  of  5  hz  is 
applied  before  convolving  the  transfer  function  with  actual  data. 

A  typical  example  of  what  happens  when  the  transfer  operators  are  applied 
to  observations  is  shown  in  Figure  9.  The  records  are  from  the  BOXCAR  event 
at  a  range  of  7.3  km.  The  vertical  and  radial  velocity  traces  are  shown  on 
the  left.  They  display  a  relatively  impulsive  and  simple  P  wave  arrival.  The 
amplitude  of  the  first  peak  is  about  30  cm. /sec.  on  either  component. 

Transfer  operators  like  those  shown  in  Figure  8  were  computed  and  convolved 
through  to  produce  the  four  derivative  traces  on  the  right.  In  this  example 
partial  derivatives  of  the  vertical  velocity  trace  are  generated  from  the 
vertical  trace  and  derivatives  of  the  radial  from  the  radial.  In  principle, 
this  need  not  be  the  case  as  we  discuss  in  the  following.  The  first  peak  of 


VELOCITY  TO  STRAIN  TRANSFER 


E-iguffi  fij.  An  example  of  the  computation  of  the  velocity  to 
strain  transfer  functions  for  the  Pahute  Mesa  structure.  The 
velocity  traces  on  the  left  are  deconvolved  from  the  strain 
traces  in  the  center  to  produce  the  transfer  operators  on  the 
right. 
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LL&U£e  2.1.  An  example  of  the  velocity  to  strain  transfer  for  one  of  the 
stations  from  the  event  BOXCAR.  The  velocity  pulses  are  shown  on  the  lv^t  and 
the  estimated  strain  pulses  on  the  right.  The  amplitude  values  refer  to  the 
first  positive  peak. 


dV/dz  is  strongly  reduced  with  respect  to  the  later  peaks.  This  is  the 
component  in  which  the  phase  shifting  by  the  strain  receiver  function  is 
generally  the  strongest.  Its  amplitude  is  also  usually  the  lowest  as  it  is  in 
this  example.  The  top  trace  is  Ezz  the  center  two  are  combined  to  produce  EZR 
and  the  bottom  is  Err  If  the  free  surface  condition  is  to  be  satisfied 
exactly,  the  center  two  traces  should  cancel.  This  clearly  will  not  occur  in 
the  example  shown.  In  order  for  this  delicate  cancellation  to  take  place,  the 
vertical  and  radial  traces  both  need  to  be  transfered  into  time  series  which 
are  exactly  proportional.  The  synthetics  are  never  really  exact,  and 
generally  the  vertical  component  of  motion  is  predicted  better  than  the 
radial.  The  radial  is  presumably  much  more  sensitive  to  lateral  velocity 
variations  directly  under  the  instrument.  In  the  following,  we  will  force  the 
free  surface  condition  to  be  satisfied  by  using  only  the  more  reliable 
vertical  records.  In  order  for  the  stress  P„  to  remain  zero 

E„--[l-2fi2/a2)E„  (8) 

Thus,  we  can  generate  one  nonzero  strain  from  the  other  in  such  a  way  that  the 
free  surface  condition  is  satisfied. 

Figure  10  shows  the  strains  along  with  the  trace  of  the  strain  tensor 
next  to  the  bottom  and  the  maximum  shear  strain  as  a  function  of  time  at  the 
bottom.  The  maximum  shear  strain  is  defined  as  the  absolute  value  of  the 
difference  between  the  largest  and  smallest  diagonal  elements  in  the 
diagonalized  strain  tensor.  The  event  BOXCAR  appears  to  have  generated  a  peak 
shear  of  74  /istrain  at  7.3  km  and  presumably  higher  at  closer  ranges.  The 
laboratory  data  in  Figure  2  extends  up  to  only  70  ^strain.  This  figure  shows 
that  if  the  waves  from  BOXCAR  propagated  through  a  material  like  sandstone, 
they  would  be  strongly  attenuated  indeed.  Figure  11  shows  the  corresponding 
stresses.  The  pressure  in  the  next  to  bottom  row  is  the  negative  of  1/3  the 


trace  of  the  stress  tensor  and  the  bottoa  trace  is  .he  maximum  shear  stress 
The  sixteen  processed  records  are  presented  in  Appendix  1  in  the  saae  format 
as  Figures  10  and  11. 

The  decay  rate  of  the  peak  shear  strain  with  range  is  shown  in  Figure  12. 
Theoretical  curves  are  shown  for  event  BOXCAR  which  is  the  largest  in  the  data 
set  and  for  SCOTCH  which  is  the  smallest.  The  curves  are  computed  assuming 
elastic  theory  is  appropriate  and  using  the  same  velocity  structure  as  in  the 
calculation  of  transfer  functions.  Note  in  Figure  7  that  this  structure  is 
consistant  with  models  derived  using  a  variety  of  different  approaches.  The 
observed  values  from  the  processed  records  are  shown  as  data  points.  The  70 
pstrain  level  is  only  half  way  up  the  vertical  scale.  The  theoretical  curve 
for  the  smallest  event  does  not  drop  below  1  ^strain  by  25  km.  Thus  the 
entire  data  set  is  within  the  strain  regime  in  which  the  nonlinear  process 
observed  in  the  laboratory  is  believed  to  be  significant.  It  is  interesting 
that  the  observations  show  the  same  rate  of  decay  as  the  theoretical  curve . 

The  curve  was  computed  using  a  theory  which  assumes  linear  elasticity  holds. 
The  nonlinear  effects,  if  they  are  indeed  significant,  should  have  dropped  the 
observed  values  below  the  elastic  curves  with  range.  However,  it  is  important 
to  remember  that  the  nonlinearity  associated  with  sliding  on  cracks  (see 
Figure  2)  probably  only  occurs  very  near  to  the  surface.  All  of  the 
generalized  rays  important  to  the  P  pulse  dive  downward  into  the  crust  and 
only  enter  the  region  where  they  might  be  attenuated  as  they  emerge  under  the 
receiver.  In  other  words,  all  of  the  signals  in  the  data  from  a  given  event 
might  be  attenuated  by  more  or  less  the  same  amount.  The  reduction  in 
amplitude  by  the  nonlinearity  would  then  be  reflected  in  an  underestimate  of 
the  absolute  size  of  the  event. 
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Figure  10,  Strains  from  a  typical  BOXCAR  station.  The  top  two 
traces  are  the  nonzero  strains  in  a  cylindrical  coordinate  system 
assuming  that  the  free  surface  condition  holds.  The  third  trace 
is  the  trace  of  the  strain  tensor,  and  the  bottom  trace  is  the 
peak  shear  strain. 


4.553xl0-1  MPa 


TIME  (SEC) 


F iqrura  11 .  Strassaa  froa  a  typical  BOXCAR  station.  The  top  two 
tracts  art  tha  nonzero  stresses  in  a  cylindrical  coordinate 
systea.  Tha  third  trace  is  pressure  and  the  fourth  is  peak  shear 
stress . 
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1 l gure  12  The  decay  of  peak  shear  strain  near  to  nuclear  explosions  The 
largest  observed  value  of  shear  strain  is  10  5  Theoretical  curves  are  showi 
for  the  largest  and  smallest  events  studied  The  observations  appeal  to 
follow  the  theoretical  predictions 


.AAA’s'-. 


STRAINS  m  STRESSES  U£AE  Ifl  EARTHQUAKES 


The  next  records  to  which  w«  will  apply  our  velocity  to  strein 
trensformetion  procedure  are  f roe  a  small  aftershock  of  the  IS  October  19 79 
Imperial  Valley  earthquake  The  aftershock  was  studied  in  some  detail  by  Liu 
and  Helmberger  (1985)  who  providad  a  mechanism,  moment  and  time  function  for 
it  They  reported  the  event  depth  as  9  5  km  Figure  13  shows  the  strong 
motion  recordings  from  it  A  map  of  the  stations,  the  aftershock  location  and 
the  mains hock  location  Is  shown  In  an  Inset  The  strong,  clear  pulse  is  the 
direct  S  wave  The  P  wave  was  only  recorded  In  its  entirety  at  a  few  stations 
and  was  too  complex  to  model  Several  polarity  changes  are  apparent  in  the 
data  Liu  and  Helmberger  (1985)  used  this  information  to  infer  that  the  even: 
had  a  vertical  strike  slip  mechanism  They  reported  a  moment  of  1  0  x  10a* 
dyne  cm  and  a  triangular  time  function  with  a  rise  of  0  1  sec  and  a  fall  of 
0  1  sec  From  our  modeling  studies,  however,  we  conclude  that  a  moment  value 
of  0  6  x  10J*  dyne  cm  and  a  source  with  a  rise  of  0  3  sec  and  a  fall  of  0  1 
sec  is  more  accurate 

Certain  unusual  characteristics  of  the  strong  motion  records  will  guide 
how  we  will  proceed  to  estimate  strain  from  velocity  in  this  Instance  Figure 
14  shows  the  three  components  of  motion  observed  at  four  of  the  stations  at  a 
representative  set  of  ranges  HOLT  is  the  closest  station  and  RRAV  the 
farthest  HOLT.  ELCE  and  BNCR  were  three  of  the  six  stations  that  recorded 
the  complete  P  as  well  as  the  S  waves  Note  in  these  records  that  the  P  wave 
is  smaller  and  much  less  coherent  than  the  S  wave  The  peak  strains  are 
carried  by  the  S  wave  pulse,  so  we  will  transform  the  S  waves  alone  and  not 
attempt  to  process  the  P  waves  Also  note  that  there  is  no  clear  SV  arrival 
on  the  vertical  records  In  some  cases  there  is  a  burst  of  incoherent  energy 
but  there  is  no  clear  long  period  pulse  as  on  the  radial  component  The  fact 
that  the  SV  pulse  is  so  small  on  the  vertical  component  can  be  explained  bv 
the  fact  that  the  shear  velocity  near  the  surface  in  Imperial  valley  is  verv 
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tisuis.  The  redial  ground  velocities  from  the  23:19,  15  October  1979 

Imperial  Valley  aftershock 


ElfeVIS  Four  typical  three  component  record  sets  from  the  Imperial  Valley 

aftershock.  The  P  wave  Is  more  complex  and  much  higher  frequency  than  the  S 
wave.  The  vertical  S  wave  component  is  small. 


low.  Thus  the  SV  ray  emerges  almost  vertically  and  the  wave  produces  little 
vertical  motion.  We  wish  to  avoid  using  the  noisy  vertical  signals  in  our 
processing  since  they  will  probably  only  degrade  the  accuracy  of  our  strain 
estimate . 

The  nonzero  strains  generated  by  an  incident  SV  wave  are  Elz  Err  and  Err 
As  before,  we  could  use  the  vertical  velocity  record  to  generate  the  first  of 
these  and  the  radial  to  generate  the  other  two.  Instead,  we  again  will  use 
relation  8  to  generate  the  information  we  need  regarding  EIt  from  Err  and  do 
it  in  such  a  way  that  the  free  surface  condition  is  automatically  satisfied. 

In  essence,  we  generate  an  estimate  of  the  vertical  signal  from  the  radial 
signal  instead  of  using  the  noisy  vertical  channel  itself.  One  other  point 
worth  noting  about  Figure  14  is  the  clear  shift  in  frequency  content  between  P 
and  S  waves.  Liu  and  Helmberger  (1985)  attribute  this  shift  to  a  low 
effective  shear  Q  in  the  Imperial  Valley.  The  value  they  used  for  shear  Q  in 
the  top  layer  of  their  crustal  model  was  only  6.2.  This  results  in  a 
relatively  distance  independent  t*  of  .132  sec.  In  the  calculations  shown  in 
the  following,  we  use  this  value  along  with  the  crustal  model  they  presented 
(see  Table  1).  A  value  for  shear  Q  as  low  as  6.2  is  surely  atypical  and  could 
easily  be  interpreted  as  an  indication  that  nonlinear  processes  like  those 
suggested  by  the  laboratory  data  shown  in  Figure  2  might  be  taking  place. 

An  example  of  the  velocity  to  strain  transfer  operators  for  the 
earthquake  case  is  given  in  Figure  15.  As  in  Figure  8,  the  relevant 
theoretical  velocity  traces  are  shown  on  the  left,  the  spatial  derivatives  of 
them  in  the  center  and  the  transfer  operators  on  the  right.  In  this  instance, 
the  velocity  traces  are  Q  (radial)  and  V  (tangential).  The  transfer  operators 
are  shown  with  the  gaussian  filter  (cutoff  5  hz)  convolved  through.  Of  the 
four  partial  derivatives  shown,  only  two  are  actually  used  in  the  calculations 
that  follow.  The  top  one  would  be  used  to  generate  Eri  but  the  other  term  in 
would  always  cancel  it  to  satisfy  the  free  surface  condition.  The  third 
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VELOCITY  TO  STRAIN  TRANSFER 

Imperial  Valley  Aftershock,  Range  =  10  km 


Flyure  15 .  An  example  of  the  computation  of  the  velocity  to  strain  transfer 
functions  for  the  Imperial  Valley  structure.  The  velocity  traces  on  the  left 
are  deconvolved  from  the  strain  traces  in  the  center  to  produce  the  transfer 
operators  on  the  right. 


transfer  function  is  always  zero  because  the  corresponding  receiver 
coefficient  is  zero  in  Table  2.  The  second  and  fourth  transfer  functions 
generate  nonzero  strain  components  Err  and  Ert.  The  range  for  the  calculations 
shown  in  Figure  15  is  10  km  just  as  in  Figure  8.  The  strains  for  the 
earthquake  are  slightly  larger  than  for  the  megaton  explosion.  This  is  not 
true  at  all  ranges  because,  as  we  shall  show  in  the  following,  the  earthquake 
strains  decay  at  a  much  slower  rate.  Also,  because  different  components  of 
the  strain  tensor  are  nonzero,  the  peak  shear  strains  are  somewhat  larger  in 
the  explosion  case.  The  transfer  operators  for  the  earthquake  source  are 
generally  simpler  than  for  the  explosion  case  meaning  that  we  are  probably 
obtaining  a  more  reliable  estimate  of  dynamic  strain  for  the  earthquake. 

Figure  16  shows  the  transfer  of  the  HOLT  velocity  record  into  strain. 

The  strain  traces  on  the  right  are  not  strongly  altered  from  the  velocity 
traces  on  the  left.  The  changes  in  polarity  that  do  or  do  not  occur  are  just 
as  predicted  in  Table  2.  The  smoothing  out  of  the  detail  is  primarily  caused 
by  the  gaussian  filter  used  in  the  deconvolution.  The  two  nonzero  strains  are 
shown  at  the  top  of  Figure  17  along  with  the  trace  of  the  strain  tensor  and 
the  maximum  shear.  It  is  of  interest  to  compare  the  peak  strains  for  the  HOLT 
record  to  those  of  the  SCOTCH  record  from  6.1  km.  The  peak  shear  for  SCOTCH 
is  about  twice  as  large  and  the  peak  compressive  strain  about  three  times  as 
large.  It  is  important  to  note,  however,  that  the  SCOTCH  source  is  actually 
much  closer  to  the  station  than  the  earthquake  source.  The  depth  of  the 
earthquake  is  9.5  km  while  that  of  SCOTCH  is  0.97  km.  A  theoretical 
calculation  of  the  strain  for  an  earthquake  source  at  the  same  depth  and  range 
as  the  SCOTCH  record  predicts  that  the  strains  from  the  earthquake  would  be 
slightly  higher.  Figure  18  shows  the  stresses  associated  with  the  HOLT 
strains.  Stresses  and  strains  from  all  16  stations  are  given  in  Appendix  2. 
Figure  19  shows  the  decay  of  peak  shear  strain  with  range.  The  observations 
are  shown  as  data  points  and  the  theoretical  prediction  of  the  model  as  a 
smooth  curve.  The  theoretical  curve  was  computed  for  the  particular 
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Fi&Vtg  lsu  An  example  of  the  velocity  to  strain  transfer  for  station  HOLT 
from  the  Imperial  Valley  aftershock.  The  velocity  pulses  are  shown  on  the 
left  and  the  estimated  strain  pulses  on  the  right. 
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strain 


Figure  19 .  The  decay  of  peak  shear  strain  near  to  the  Imperial  Valley 
aftershock . 
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azimuth  of  45'  Here  the  SH  and  SV  waves  are  equal  amplitude  and  the  peak 
shear  strain  is  relatively  high  The  data  come  from  a  variety  of  azimuths 


which  is  why  many  of  the  data  points  fall  under  the  curve  Also,  the  sout^e 
model  of  Liu  and  Helmberger  (1985)  fails  to  correctly  predict  the  ratio  o{  s. 
to  SH.  The  SH  data  alone  suggest  a  moment  of  42  x  102‘  dyne  cm  and  the  S'. 

.73  x  102*  dyne-cm.  The  SH  wave  apparently  has  a  more  important  effect  in 
determining  the  peak  shear  strain  in  many  cases  The  peak  shear  for  the 
earthquake  at  the  surface  of  the  earth  is  lower  than  for  the  explosions  but 
decays  much  less  rapidly  with  range  This  is  because  of  the  vertical 
radiation  pattern  of  the  source  and  because  the  earthquake  source  is  deeper 
The  strains  generated  by  the  earthquake  are  large  enough  so  that  the  nonlinear 
effects  illustrated  in  Figure  2  are  potentially  significant 

DISCUSSION 

The  bulk  of  all  seismic  observations  have  been  very  successfuly  explained 
using  linear  elastic  theory.  It  seems  doubtful  that  nonlinear  processes  could 
be  of  great  importance  without  having  been  noted  previously  On  the  other 
hand,  the  science  is  still  evolving  and  new  types  of  data  are  being  studied 
Very  high  frequencies  (5  to  20  hz)  are  being  studied  for  potential  use  in 
discrimination  between  earthquakes  and  explosions.  A  thin,  shallow  layer  with 
nonlinear  response  might  have  an  effect  on  the  generation  of  such  energy  but 
would  be  of  no  significance  to  1  hz  or  lower  frequency  energy  If  the 
nonlinear  zone  has  different  characteristics  for  explosions  and  earthquakes, 
it  might  alter  how  regional  phases  are  initiated  in  the  two  cases  Such  a 
layer  would  obviously  also  be  important  to  the  generation  of  free  surface 
phases  such  as  pP.  At  low  frequency  such  phases  would  appear  as  elastic 
reflections,  whereas  at  very  high  frequency  they  would  appear  to  be  strongly 
attenuated.  Some  observations  of  nuclear  explosions  suggest  that  this  is  the 


case  Another  piece  of  evidence  that  toae  nonlinear  losses  are  occurring  is 
the  strong  motion  data  of  Liu  and  Helaberger  (1985)  which  does  suggest  that 
effective  Q  in  the  near  field  of  che  Imperial  Valley  earthquake  was  very  low 


It  is  Important  to  note  that  even  if  the  nonlinear  process  indicated  by 
the  laboratory  data  does  occur  in  the  earth  and  is  significant,  it  is  not 
clear  exactly  how  it  would  manifest  itself  At  one  level  of  approximation 
effective  Q  could  be  considered  to  be  a  function  of  time  with  a  value  dictated 
by  the  strain  wave  field  It  is  difficult  to  imagine  exactly  what  effect  this 
would  have  Stewart  et  al  (1983)  have  suggested  a  model  for  how  the 
nonlinearity  would  depend  on  the  density  of  cracks  in  the  medium  and  the 
overburden  pressure  presuming  that  it  is  indeed  related  to  frictional  sliding 
on  cracks  They  proposed  that 


Q-P*''/kkE  (9) 

where  P  is  the  overburden  pressure,  k  is  a  constant  function  of  the 
material's  elastic  parameters,  X  is  the  crack  density  and  E  is  the  strain 
amplitude  Day  and  Minster  (1986)  suggest  an  equivalent  linear  method  for 
solving  wave  propagation  problems  in  materials  behaving  in  a  weakly  nonlinear 
fashion,  but  their  method  does  not  adapt  easily  to  realistic  media.  Much 
additional  progress  will  be  needed  before  the  role  played  by  high  strain 
nonlinearity  in  seismic  wave  propagation  is  understood. 

CONCLUSIONS 

The  close  relationship  between  velocity  and  strain  wave  pulses  in  a  whole 
space  appears  to  be  maintained  for  the  most  part  in  a  layered  half  space.  The 
large  data  bases  of  near  field  velocity  records  which  have  been  collected  over 
the  years  can  thus  be  transformed  into  a  data  base  of  near  field  dynamic 
stress  and  strain  records.  Near  field  strains  for  a  large  explosion  appear  to 
be  as  high  as  10‘3  at  the  surface  of  the  earth.  Those  near  to  a  small 
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